\documentclass[a4paper,12pt,titlepage]{ctexart}
\usepackage{xeCJK}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{float}
\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{fancyhdr}
\usepackage{caption}
\usepackage{textcomp}
\usepackage[left=2.5cm,right=2.5cm,top=2.5cm,right=2.5cm]{geometry}

\pagestyle{headings}
\title{程序设计报告}
\author{郭畅}
\begin{document}
    \maketitle
\section{问题}
求解方程
\begin{equation}\label{eq:1}
\frac{\partial^2}{\partial s^2}\left(\frac{\xi(s)}{h(s)}\right)+\frac{\partial}{\partial s}\ln\left\{ h(s)^2B(s) \right\}\frac{\partial}{\partial s}\left(\frac{\xi(s)}{h(s)}\right)+\left(2\pi f\frac{\sqrt{\mu_0\rho}}{B(s)}\right)^2\left(\frac{\xi(s)}{h(s)}\right)=0
\end{equation}
其中有：
\begin{align*}
&h(s)=r\sin\theta\\
&B(L,\theta)=\frac{B_0}{L^3}\frac{\sqrt{1+3\cos^2\theta}}{\sin^6\theta}\\
&s(L,\chi=\cos\theta)=\frac{LR_E}{2}\left[\chi\sqrt{1+3\chi^2}+\frac{1}{\sqrt{3}}\ln\left(\sqrt{1+3\chi^2}+\sqrt{3}\chi\right)\right]
\end{align*}
磁力线表达式为\[ r=LR_E\sin^2\theta \]
边界条件为\begin{align*}
&\xi(s_{m1})=\xi(s_{m2})=0\\
&s_{m1}\text{和}s_{m2}\text{由}r=LR_E\sin^2\theta\equiv R_E\text{确定}
\end{align*}
取 $ L=6.6,f=0.019,0.051,0.082,0.113,0.144,0.175 $\\
使用到的常数为
\begin{align*}
&R_E=6376\times 10^3 \mathrm{m}\\
&B_0=3.12\times 10^{-5} \mathrm{T}\\
&\mu_0=4\pi\times 10^{-7}\mathrm{H}/\mathrm{m}\\
&\rho=1\times 10^{-6}\times 1.6726219\times 10^{-27} \mathrm{kg}/\mathrm{m}^3
\end{align*}
\section{算法分析}
\subsection{方程}
我们令\[\left\lbrace \begin{aligned}
y(s)&=\frac{\xi(s)}{h(s)}\\
p(s)&=\frac{\partial}{\partial s}\ln\left\{ h(s)^2B(s) \right\}\\
q(s)&=\left(2\pi f\frac{\sqrt{\mu_0\rho}}{B(s)}\right)^2
\end{aligned}\right. \]
得到简化的方程
\begin{equation}\label{eq:2}
\frac{\partial^2 y(s)}{\partial s^2}+p(s)\frac{\partial y(s)}{\partial s}+q(s)y(s)=0
\end{equation}
由于$ h(s),B(s) $以$ \theta $为自变量表示比较方便，因此将方程变换为以$ \theta $为自变量。
\begin{align*}
\frac{\mathrm{d}y}{\mathrm{d}s}&=\frac{\mathrm{d}y}{\mathrm{d}\theta}\frac{\mathrm{d}\theta}{\mathrm{d}s}\\
&=\frac{\mathrm{d}y}{\mathrm{d}\theta}/\frac{\mathrm{d}s}{\mathrm{d}\theta}\\
&=\frac{1}{s'_\theta}\frac{\mathrm{d}y}{\mathrm{d}\theta}\\
\frac{\mathrm{d}^2y}{\mathrm{d}s^2}&=\frac{\mathrm{d}}{\mathrm{d}s}\left(\frac{1}{s'_\theta}\frac{\mathrm{d}y}{\mathrm{d}\theta}\right)\\
&=\frac{1}{s'_\theta}\left(\frac{1}{s'_\theta}\frac{\mathrm{d}y}{\mathrm{d}\theta}\right)\\
&=\frac{1}{s'_\theta}\frac{y''_\theta s'_\theta-y'_\theta s''_\theta}{(s'_\theta)^2}\\
&=\frac{y''_\theta s'_\theta-y'_\theta s''_\theta}{(s'_\theta)^3}
\end{align*}
代入方程(\ref{eq:2})有
\begin{gather*}
\frac{1}{(s'_\theta)^2}y''_\theta+\left(\frac{p(s(\theta))}{s'_\theta}-\frac{s''_\theta}{(s'_\theta)^3}\right)y'_\theta+q(s(\theta))y(\theta)=0\\
y''_\theta+\left(p(s(\theta))s'_\theta-\frac{s''_\theta}{s'_\theta}\right)y'_\theta+q(s(\theta))(s'_\theta)^2y(\theta)=0
\end{gather*}
令\[\left\{\begin{aligned}
P(\theta)&=p(s(\theta))s'_\theta-\frac{s''_\theta}{s'_\theta}\\
Q(\theta)&=q(s(\theta))(s'_\theta)^2
\end{aligned}\right.\]
得以$ \theta $为自变量的方程
\begin{equation}\label{eq:3}
y''_\theta+P(\theta)y'_\theta+Q(\theta)y(\theta)=0
\end{equation}
令\[z(\theta)=y'_\theta\]
得方程组
\begin{equation}\label{eq:4}
\left\lbrace\begin{aligned}
&\frac{\mathrm{d}y}{\mathrm{d}\theta}=z\\
&\frac{\mathrm{d}z}{\mathrm{d}\theta}=-Pz-Qy
\end{aligned}\right.
\end{equation}
这可以采用龙格-库塔法求解。
\subsection{边界条件}
求解$ LR_E\sin^2\theta=R_E $，得$ \sin\theta=\pm\frac{1}{\sqrt{L}} $\\
磁力线分布如图
\begin{figure}[H]
    \centering
    \includegraphics[width=0.6\linewidth]{mag_line.eps}
    \caption{磁力线分布}
\end{figure}
图中蓝色为磁力线，红色为地球表面。易知只需求$ \sin\theta=\frac{1}{\sqrt{L}} $在$ (0,\pi) $上的两个解即可。\\
设解为$ \theta_1,\theta_2 $，相应的s使$ \xi(s)=0 $，从而\[
y(\theta_{1,2})=\frac{\xi(s)}{h(s)}=\frac{\xi}{r(\theta_{1,2})\sin\theta_{1,2}}=\frac{\sqrt{L}\xi}{R_E}=0 \]
\section{算法设计}
\subsection{大概步骤}
按照以上思路，可以整理出算法的大概步骤
\begin{enumerate}
    \item 对给定的\textit{L}，求解$\theta$的范围$ (\theta_1,\theta_2) $，d得到边界条件为$ y(\theta_1)=y(\theta_2)=0 $
    \item 划分网格，生成自变量的数组Theta
    \item 按照方程组(\ref{eq:4})求解方程，得到$ y(\theta) $
    \item 计算$ h(\theta)=h(s(\theta))=LR_E\sin^3\theta $
    \item 求出$ \xi(\theta)=y(\theta)h(\theta) $
    \item 求$ s=s(\theta) $
    \item 保存$ \xi\sim s $的数据
\end{enumerate}
\subsection{一些细节}
\subsubsection{$ s'_\theta $和$ s''_\theta $}
因为\[s(\chi=\cos\theta)=\frac{LR_E}{2}\left[\chi\sqrt{1+3\chi^2}+\frac{1}{\sqrt{3}}\ln\left(\sqrt{1+3\chi^2}+\sqrt{3}\chi\right)\right]\]
所以有
\[ s'_\theta=s'_\chi\chi'_\theta=-s'_\chi\sin\theta \]
记$ \beta=\sqrt{1+3\chi^2} $，则$ \beta'_\chi=3\chi/\beta $
\begin{align*}
s'_\chi&=\left\{\frac{LR_E}{2}\left[\chi\beta+\frac{1}{\sqrt{3}}\ln(\beta+\sqrt{3}\chi)\right]\right\}'_\chi\\
&=\frac{LR_E}{2}\left[\beta+\frac{3\chi^2}{\beta}+\frac{1}{\sqrt{3}}\frac{\frac{3\chi}{\beta}+\sqrt{3}}{\beta+\sqrt{3}\chi}\right]\\
&=\frac{LR_E}{2}\left[\beta+\frac{\beta^2-1}{\beta}+\frac{1+\frac{\sqrt{3}\chi}{\beta}}{\beta+\sqrt{3}\chi}\right]\\
&=\frac{LR_E}{2}\left[\beta+\beta-\frac{1}{\beta}+\frac{1}{\beta}\right]\\
&=LR_E\beta\\
s''_\chi&=\frac{3LR_E\chi}{\beta}
\end{align*}
即
\begin{gather}
\label{ds} s'_\chi=LR_E\sqrt{1+3\chi^2} \\ 
\label{dds} s''_\chi=\frac{3LR_E\chi}{\sqrt{1+3\chi^2}} 
\end{gather}
二阶导数
\begin{align*}
s''_\theta&=(s'_\theta)'_\theta\\
&=(-\sin\theta s'_\chi)'_\theta\\
&=-s'_\chi\cos\theta-\sin\theta s''_{\chi\theta}\\
&=-s'_\chi\cos\theta+\sin^2\theta s''_\chi
\end{align*}
因此
\begin{align*}
\frac{s''_\theta}{s'_\theta}&=\frac{-s'_\chi\cos\theta+\sin^2\theta s''_\chi}{-s'_\chi\sin\theta}\\
&=\frac{\cos\theta}{\sin\theta}-\frac{s''_\chi}{s'_\chi}\sin\theta\\
&=\frac{\cos\theta}{\sin\theta}-\frac{3\chi}{\beta^2}\sin\theta\\
&=\frac{\cos\theta}{\sin\theta}-\frac{3\chi}{1+3\chi^2}\sin\theta\\
&=\frac{\cos\theta}{\sin\theta}-\frac{3\cos\theta\sin\theta}{1+3\cos^2\theta}
\end{align*}
\subsubsection{$ P(\theta) $}
已知$ h(s)=r\sin\theta=LR_E\sin^3\theta $和$ B(\theta)=\frac{B_0}{L^3}\frac{\sqrt{1+3\cos^2\theta}}{\sin^6\theta} $\\
有
\begin{align*}
p(s)&=\frac{\partial}{\partial s}\ln\left\{ h(s)^2B(s) \right\}\\
&=\frac{\partial}{\partial s}\left(2\ln h(s)+\ln B(s)\right)\\
&=\frac{1}{s'_\theta}\frac{\partial}{\partial\theta}\left(2\ln h+\ln B\right)\\
&=\frac{1}{s'_\theta}\left(2\frac{h'_\theta}{h}+\frac{B'_\theta}{B}\right)\\
&=\frac{1}{s'_\theta}\left(\frac{6\cos\theta}{\sin\theta}-\frac{3\sin\theta\cos\theta}{1+3\cos^2\theta}-\frac{6\cos\theta}{\sin\theta}\right)\\
&=-\frac{3\sin\theta\cos\theta}{s'_\theta(1+3\cos^2\theta)}
\end{align*}
故\begin{align*}
P(\theta)&=p(s(\theta))s'_\theta-\frac{s''_\theta}{s'_\theta}\\
&=-\frac{3\sin\theta\cos\theta}{1+3\cos^2\theta}-\left(\frac{\cos\theta}{\sin\theta}-\frac{3\cos\theta\sin\theta}{1+3\cos^2\theta}\right)\\
&=-\frac{\cos\theta}{\sin\theta}\\
&=-\cot\theta
\end{align*}
\subsubsection{$ Q(\theta) $}
已知\[ q(s)=\left(2\pi f\frac{\sqrt{\mu_0\rho}}{B(s)}\right)^2 \]
因此\begin{align*}
Q(\theta)&=q(\theta)(s'_\theta)^2\\
&=\left(2\pi f\frac{\sqrt{\mu_0\rho}}{B}s'_\theta\right)^2\\
&=\left(2\pi f\frac{L^3\sqrt{\mu_0\rho}}{B_0}\frac{\sin^6\theta}{\sqrt{1+3\cos^2\theta}}LR_E\beta\sin\theta\right)^2\\
&=\left(2\pi f\frac{L^4R_E\sqrt{\mu_0\rho}}{B_0}\sin^7\theta\right)^2\\
&=Q_0\sin^{14}\theta
\end{align*}
方程(\ref{eq:3})可写为：$ y''-y'\cot\theta+Q_0y\sin^{14}\theta=0 $
\end{document}
